Optimal Portfolio Allocation

An investment universe of the following risky assets with a dependence structure (correlation) is given:

\begin{equation*} \begin{matrix} \textbf{Asset} & \boldsymbol\mu & \boldsymbol\sigma & \boldsymbol\omega\\ A & 0.02 & 0.05 & \omega_1\\ B & 0.07 & 0.12 & \omega_2\\ C & 0.15 & 0.17 & \omega_3\\ D & 0.20 & 0.25 & \omega_4 \end{matrix} \qquad \qquad R = \begin{pmatrix} 1 & 0.3 & 0.3 & 0.3\\ 0.3 & 1 & 0.6 & 0.6\\ 0.3 & 0.6 & 1 & 0.6\\ 0.3 & 0.6 & 0.6 & 1 \end{pmatrix} \end{equation*}

Question 1. Consider minimum variance portfolio with a target return m.

\begin{equation*} \underset{\boldsymbol\omega}{\arg\min} \frac{1}{2}\omega'\Sigma\omega \qquad s.t. \quad \omega'1 = 1 \quad \mu_\Pi = \omega' \mu = m \end{equation*}

• Formulate the Lagrangian and give its partial derivatives.
• Write down the analytical solution for optimal allocations $ \omega^* $ (derivation not required).
• Inverse optimisation: generate above 700 random allocation sets (vectors) 4 × 1, those will not be optimal allocations.
  • Standardise each set to satisfy $\omega'1 = 1$, in fact you can generate 3 allocations and compute the 4th.
  • For each vector of allocations compute $\mu_\Pi = \omega'\mu$ and $\sigma_\Pi = \sqrt{\omega'\Sigma\omega}$.
  • Plot the cloud of points of $\mu_\Pi$ vertically on $\sigma_\Pi$ horizontally. Explain this plot.

Hint: Since we treat this as an inverse optimisation, there is no computation of $\omega^*$ from the ready formula.

Answers and code

Form the Langrange function with two multipliers $\lambda$ and $\gamma$ for a given vector of weights $w$ \begin{equation*} L(w, \lambda, \gamma) = \frac{1}{2}w'\Sigma w + \lambda(m - w'\mu) + \gamma(1 - w'1) \end{equation*}

where $(m - w'\mu)$ is the return constraint and $(1 - w'1)$ is the budget constraint.

First order condition

The optimal portfolio is the point where the gradient of $L(w, \lambda, \gamma) = 0$.

Therefore $\frac{\partial L}{\partial w} = 0$.

\begin{equation*} \frac{\partial L}{\partial w} = \frac{1}{2} 2 \Sigma w + \lambda (-\mu) + \gamma(-1) = \Sigma w - \lambda \mu - \gamma 1 = 0 \end{equation*}

Second order condition

Given this is a minimization problem, we seek the Hessian to be greater than 1

\begin{equation*} \frac{\partial^2 L}{\partial w^2} = \Sigma \end{equation*}

As the covariance matrix is positive, the 2nd order condition is satisfied.

Therefore the optimal allocation $w^*$ can be derived

\begin{equation*} \Sigma w - \lambda \mu - \gamma 1 = 0 \end{equation*}\begin{equation*} w^* = \Sigma^{-1} (\lambda \mu + \gamma 1) \end{equation*}

Substituting the above optimal weight $w^*$ into the constraint

\begin{equation*} \mu'w = \mu' \Sigma^{-1} (\lambda \mu + \gamma 1) = \mu' \Sigma^{-1} \lambda \mu + \mu' \Sigma^{-1}\gamma 1 = m \end{equation*}

and \begin{equation*} 1'w = 1' \Sigma^{-1} (\lambda \mu + \gamma 1) = 1' \Sigma^{-1} \lambda \mu + 1' \Sigma^{-1}\gamma 1 = 1 \end{equation*}

Setting scalars for ease of presentation and computation \begin{equation*} A = 1\Sigma^{-1} 1' \end{equation*}

\begin{equation*} B = \mu' \Sigma^{-1} 1 = 1' \Sigma^{-1} \mu \end{equation*}\begin{equation*} C = \mu' \Sigma^{-1} \mu \end{equation*}

The constraint functions then become

\begin{equation*} C \lambda + B \gamma = m \end{equation*}\begin{equation*} B \lambda + A \gamma = 1 \end{equation*}

This is a system of 2 equations with 2 unknowns

\begin{equation*} \lambda = \frac{m - B \gamma}{C} \end{equation*}

therefore

\begin{equation*} \frac{Bm - B^2 \gamma}{C} + A \gamma = 1 \end{equation*}\begin{equation*} Bm - B^2 \gamma + A C \gamma = C \end{equation*}\begin{equation*} A C \gamma - B^2 \gamma = C - Bm \end{equation*}\begin{equation*} \gamma = \frac{C - Bm}{AC - B^2} \end{equation*}

in the same manner

\begin{equation*} \lambda = \frac{Am - B}{AC - B^2} \end{equation*}

Finally, substituting back into optimal weight $w^*$

\begin{equation*} w^* = \Sigma^{-1} (\lambda \mu + \gamma 1) = \frac{1}{AC - B^2} \Sigma^{-1}[(A \mu - B1)m + (C1 - B \mu)] \end{equation*}

# Import the relevant packages
# Import numpy
import plotly.express as px
import numpy as np
from numpy import *
from numpy.linalg import multi_dot

# Import pandas
import pandas as pd

# Import plotly express
px.defaults.template = 'ggplot2'
px.defaults.width, px.defaults.height = 600, 400
import plotly.io as pio
pio.renderers.default='notebook'

# Set the matrices for mu and sigma
mu = np.array([0.02, 0.07, 0.15, 0.2]).reshape((4, 1))
sigma = np.array([0.05, 0.12, 0.17, 0.25]).reshape((4, 1))
R = np.array([[1, 0.3, 0.3, 0.3], [0.3, 1, 0.6, 0.6], [
             0.3, 0.6, 1, 0.6], [0.3, 0.6, 0.6, 1]]).reshape(4, 4)

# Compute the covariance matrix usings SRS where S = the diagonal standard deviation matrix
S = np.diag(sigma.flatten())
covariance_matrix = multi_dot([S, R, S])

# Set values for N and initiate an empty list for the
N = 10000
weights_arr = []

# Generate N random weights for the 4 assets
for i in range(N):
    weights = random.random(4)
    # Standardise these weights
    weights /= sum(weights)
    weights_arr.append(weights)

# Add them to a dataframe
df = pd.DataFrame(weights_arr, columns=['w_1', 'w_2', 'w_3', 'w_4'])

df.head()

	w_1	w_2	w_3	w_4
0	0.229252	0.159475	0.305795	0.305478
1	0.016505	0.413196	0.182878	0.387420
2	0.085283	0.621982	0.237425	0.055310
3	0.204154	0.141498	0.206525	0.447823
4	0.325325	0.192920	0.286947	0.194808

df.tail()

	w_1	w_2	w_3	w_4
9995	0.270300	0.261418	0.429156	0.039126
9996	0.251775	0.355114	0.029828	0.363283
9997	0.141022	0.059185	0.587882	0.211911
9998	0.382156	0.239544	0.334225	0.044075
9999	0.623285	0.091566	0.218903	0.066245

# loop through weights, extract matrix, calculate portfolio mean and sd, append back to df
for index, row in df.iterrows():
    # convert to a matrix of weights
    w = np.array(row).reshape(1, 4)
    # Compute mean for the portfolio, add it to the dataframe
    mu_pi = multi_dot([w, mu])
    df.loc[index, 'mu_pi'] = mu_pi.flatten()
    # Compute the sd for the portfolio, add it to the dataframe
    sigma_pi = np.sqrt(multi_dot([w, covariance_matrix, w.T]))
    df.loc[index, 'sigma_pi'] = sigma_pi.flatten()

df.head()

	w_1	w_2	w_3	w_4	mu_pi	sigma_pi
0	0.229252	0.159475	0.305795	0.305478	0.122713	0.133261
1	0.016505	0.413196	0.182878	0.387420	0.134170	0.155308
2	0.085283	0.621982	0.237425	0.055310	0.091920	0.115161
3	0.204154	0.141498	0.206525	0.447823	0.134531	0.151253
4	0.325325	0.192920	0.286947	0.194808	0.102015	0.110853

df.tail()

	w_1	w_2	w_3	w_4	mu_pi	sigma_pi
9995	0.270300	0.261418	0.429156	0.039126	0.095904	0.107151
9996	0.251775	0.355114	0.029828	0.363283	0.107024	0.129425
9997	0.141022	0.059185	0.587882	0.211911	0.137528	0.145700
9998	0.382156	0.239544	0.334225	0.044075	0.083360	0.093486
9999	0.623285	0.091566	0.218903	0.066245	0.064960	0.073557

# Plotting the results
fig = px.scatter(df, x='sigma_pi', y='mu_pi', color='mu_pi', color_continuous_scale='purp',
                 labels={'sigma_pi': 'Expected Volatility',
                         'mu_pi': 'Expected Return'},
                 title='Efficient Frontier highlighting maximum and minimum return')
idx_max_mu = df['mu_pi'].idxmax()
idx_min_mu = df['mu_pi'].idxmin()

fig.add_scatter(x=[df['sigma_pi'].loc[idx_max_mu]], y=[df['mu_pi'].loc[idx_max_mu]], marker=dict(color='green', size=10, symbol='cross'),
                name='Maximum Return').update(layout_showlegend=False)
fig.add_scatter(x=[df['sigma_pi'].loc[idx_min_mu]], y=[df['mu_pi'].loc[idx_min_mu]], marker=dict(color='red', size=10, symbol='cross'),
                name='Minimum Return').update(layout_showlegend=False)
fig.show()

The above plot builds up an example through simulated data showing the efficient frontier of a portfolio of 4 risky assets.